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ABSTRACT 


During hot firing ground testing of the Space Shuttle's 
Main Engine, vibrations of the liquid oxygen pump have 
occurred at frequencies which cannot be explained by the 
linear Jeffcott model of the rotor. The model becomes 
nonlinear after accounting for deadband, side forces, 
and rubbing. Two phenomena present in the numerical 
solutions of the differential equations are unexpected 
periodic orbits of the rotor and "tracking" of the 
nonlinear frequency. A multiple-scale asymptotic 
expansion of the differential equations is used to give 
an analytic explanation of these characteristics. 
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1. Introduction: 


Excessive vibrations of the liquid oxygen turbopump 
in the Space Shuttle's Main Engine have been recorded 
during hot firing ground testing. Examination of the 
power spectral density (PSD) plot reveal frequencies 
which cannot be explained using the linear rotordynamics 
model of Jeffcott [5]. Consequently, numerous investi- 
gations have been undertaken to study such rotors and to 
provide descriptions of the solution of the two, coupled, 
second-order differential equations which describe the 
motion of the rotor's center of mass. Following the early 
work in rotordynamics by Yomamoto [7] , one introduces a 
nonlinearity to the Jeffcott equations by including the 
effect of bearing clearance or deadband. In the pump this 
deadband refers to the load carriers (ball bearings) and 
physically describes the clearance between the outer race 
of the bearing and the support housing. Both empirical 
results by Childs [1,7] and Gupta et al . [4] and numerical 
solutions using the fourth order Runge-Kutta algorithm by 
Control Dynamics Company [3] have been helpful in under- 
standing the rotor's motion for the nonlinear problem. 

This report extends the earlier work by using analytic 
expressions obtained from singular asymptotic expansions 
(methods of multiple scales) to quantize the solution. 

Section 2 presets the general analysis of a Jeffcott 
rotor with deadband and sinusoidal forcing. This discus- 
sion includes a derivation of the characteristic, non- 
linear frequency which is responsible for "tracking", a 
shift in the subsynchronous frequency as a function of the 
external force. Then the method of multiple scales provides 
justification for inclusion or exclusion of the subsynchro- 
nous term in the solution expression. 

Section 3 displays three examples. The first example 
illustrates the major theoretical devices of section 2. 
Example 2 uses data from CDC report [3] where the subsyn- 
chronous frequency begins at one-half the forcing frequency. 
It reiterates the mathematical development of section 2. 

The last example uses the data of example 1 to produce a 
frequency- response plot for the nonlinear problem. How 
this plot differs from the corresponding linear case and 
what one may expect in the general case are the major points 
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made in this example. All the examples use the method of 
Runge-Kutta to obtain numerical solutions. 

Section 4 concludes the report with brief descriptions 
of extension and related problems. These include multiple 
external forces (e.g. sideforce) , rubbing, asymmetric stiff- 
ness, and stability analysis. 

2 . General Theory : 

The linear Jeffcott equations which describe the motion 
of a rotor in the inertial, Cartesian coordinate system are 
these : 

(1.) my=-Csy-K s y-QsZ+muw : ’coswt 

(2.) mz=-Csz+Q s y -K s z+muw 2 sinu)t 
where 


m - mass (kg.) 

- seal damping (kg./s.) 

K s = seal stiffness (kg./s. z ) 

Qs= cross-coupling stiffness of seal (kg./ s. 2 ) 

u = displacement of the shaft center of mass 
from the geometric center (m.) 

w = angular velocity of the shaft (rad./s. ) 

For the model to include bearing forces which hold the 
rotor in position, one adds tie terms 

-KsCy-yfi/ /y z + z 2 ) + mKb ( z - z 6 / /y ^ + z 2 ) 


and 

-uK B (y-yS//y 7 +P') -Kg (z-z6//y 2 +z 2 ) , 
respectively, to the right-hand sides of equations (1.) 
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and ( 2 .)* Here Kb (kg./s. 2 ) is the bearing stiffness, 6 
(m.) is the clearance or deadband between the housing and 
the bearing, and y (nondimensional) is the coefficient of 
friction between the housing and the beaiing. These bear- 
ing forces occur only when /y 2 + z 2 > 6 ; otherwise, they are 
zero. Equations (l.)-(2.) then become 

(3.) y+ (C s /m)y+ ( 1/m) [K^+KgCl-fi/r) ]y+(l/m) [Q s -yK B (l -6/r) ] z 
=uw 2 cosa)t 

(4 . ) z+(C s /m) z- (1/m) [Q s -uK B (l-6/r) ]y+(l/m) [K S +K B (1 - 6/r) ] z 
=U(o 2 sinwt 

when r = /y 2 + z > 6 ; otherwise, K B =0 . 

Equations (3 )-(4.) can be put in nondimensional form using 
a displacement g and a frequency a. Thus, using Y=y/g, 
Z=z/g, and r=at, the dimensionless equations are these: 

(5.) Y ' ’+CY '+ [A+k(l -A/R] Y+ [B-yk (1 -A/R) ] Z = E4> 2 co.-<n 

(6.) Z ' ’+CZ [B-yk(l-A/R) ]Y+[A+k(l-A/R) ]Z = E4> 2 sin4>x 

where prime denotes differentiation with respect to r and 
C=C s /m/o , A=K s /m/o 2 , k=K B /m/o 2 , B=Q s /m/a 2 , A=6/g , R=r/g, 
E=u/g, and 4 >=w/a. 

Equations (5.) -(6.) can be reduced to the following single 
equation by defining W=Y+iZ: 

(7.) W '+CW '+{A+k(l - A/ | W | ) + i [ - B+yk (1 * A / | W| ) ] }W=E<J> 2 exp(i4>T) 

Furthermore, the polar form of equations (5.) -(6.) is 

(8.) R ’ '+CR , + [A+k(l-A/R)-(O , ) 2 ]R=E<}) 2 cos((})T-0) 

(9.) R0 ' ' + (2R'+CR)0' = R[B-yk(l-A/R)]+E4) 2 sin(4)T-0) 

where R=(Y 2 +Z 2 ) ^ 2 and 0=Arctan (Z/Y) . Since y is nondimen- 
sional and typically small, one may regard it as zero with- 
out affecting the qualitative results. 

Consider the homogeneous (E=0) equation corresponding 
to equation (7.). If this equation were also linear (A=0) , 


then exponentially growing or decaying solutions would 
result for a given set of system parameters In the special 
case that (B/C) 2 *A+k, a sinusoidal solution is obtained with 
frequency 6=B/C. To see this, consider the characteristic 
equation for W=exp(mi): 

m 2 +Cm+ [A+k-iB] =0 

m=-C/2±{C 2 /4-A-k + iB} l/2 

m=-C-iB/C, iB/C. 

In the nonlinear, homogeneous problem, k is replaced by 
k(l-A/R); hence, if R is a constant, then there is a wide 
spectrum for which (B/C) 2 =A+k(l-A/R) ; i.e., if 

(10.) A< (B/ C) 2 <A+k , 

then there is a constant value of R (with R>A) for which 

(B/C) 2 =A+k (1 -A/R) . This value of R is denoted by a and 

the corresponding frequency by 6o=B/C. This frequency is 

labeled the characteristic, nonlinear frequency. Thus, 

whenever inequality (10.) is satisfied, equations (5.)-(6.) 

with E=0 have steady-state solutions Y=a cos(Bot) and Z = a sin (B 0 x) . 


The same results can be found from the polar equations 
(8.) -(9.) by assuming R and O' are constants. Then with E=0 
those equations become 


[A+k(l-A/R)-(0 *) 2 ] R=0 , C 0 ' =B • 

Thus, 8 o =0 ' = B/C and a = R=kA/ (A+k- 6 o 2 ) . 

Notice that 8 o =B/C£A+k = { (K s +Kb) / m} /o=wo , the dimen- 

sionless natural frequency of the linear system. Thus, in 
considering the general nonhomogeneous problem, it is nec- 
essary to be aware of these three dimensionless frequencies: 


Bo 

O ) 0 


the characteristic, nonlinear frequency, 
the natural frequency, 
the driving frequency. 


% 
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Either So or u) 0 is an appropriate choice for o, the non- 
dimensionalizing frequency. Correspondingly, one would 
select either a (with So)or u (with w 0 ) as the base dis- 
placement g. 

One final rearrangement of equation (9.) is made here 
to emphasis the characteristic, nonlinear frequency: 

(11.) W* , +CW' + <W=ef (W)+E4> 2 exp(i<J>x) 

where <~A+k(l-A/a) +i \ -B+yk(l -A/a) ] and 

f(W)*kA(l+iy)/e[l/|W|-l/a]W. 

Asymptotic expansions may now be introduced to solve 
equation (11.). The singular method of multiple scales, 
as described by Nayfeh [6], is appropriate for this problem. 
The method of averaging gives similar results. 

Instead of one time scale x, assume the problem depends 
on many time scales: 

To=T. Ti=ET, T2=£ 2 T,... 

Henceforth, only T 0 and Ti are used. Let W(x) =W(T 0 ,Ti) = 

Wo (To ,Ti) +eWi (To , Ti)+... Equation (11.) becomes a partial 
differential equation since 

d/di=(3/8T 0 ) (dTo/dr)+(3/3Ti) (dTi/dx) =D 0 +eDi 

and (d 2 /dx 2 ) =(Do+eDi) 2 . 

Thus, one finds 

(12.) (D 0 + cD 1 ) 2 (^o+£W 1 + . . .)+C(Do+eDi) (W 0 +eW! + . . .) 

+ <(W 0 + eWi+. . ,)=ef(Wo+eWi+. . .)+E<f> 2 exp(i<f>To) • 

Equating like powers of e yields 

(13 . ) e 0 : Do 2 W 0 +CDo Wo + <Wo = E0 2 exp( i0T 0 ) • 

This is a linear problem with this steady-staxe solution 
W 0 =Mexp ( IB oTo ) + Nexp (i<f>T 0 ) 
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N 


i 

where N=*E4> 2 / (-<j) 2 + iC4» + <) and M=*M(Ti). To determine M one 
must examine the e-order problem and chose M to eliminate 
secular terms; see Nayfeh [6]: 

e 1 : D 0 1 W 1 +CWi+kW 1 != -2DoDiWo -CDiW 0 + f (W 0 ) . 

with L*kA/e, the right-hand side of the last equation 
becomes 

-2 iBoM T exp ( iB 0 T 0 ) -CM' exp ( iB 0 T 0 ; 

+ L ( 1 / | W o | -1/a) [Mexp(iBoTo) +Nexp ( i <?T o ) ] 
where |W 0 |={|M| 2 +|N| 2 +RNexp[ i(<J>-6o)T 0 ]+MNexp[ i(Bo -4>)To ] } 1/2 

To avoid secular terms one requires that the collective co- 
efficient of exp (iBoTo) be zero. Although an analytic 
solution of the differential equation for M(Ti) has not 
been found, one can qualitatively assess M based on a 
similar problem (van der Pol's equation) and specific 
numerical results (presented in the next section) . 

Since M(Ti) is complex, it may be written as 
M(Ti)=p(Ti)exp[i&(Tx)] . Thus, 

Wt = p (Ti) exp [ iB oT : : iB (Tx) ] +Nexp (i4>T o) or, assuming 
B(Ti) is analytic near t = 0, W 0 = p (T i) exp [ i (B o + eB i) t + . . . ] | + 

Nexp(i<}>T). Thus the fundamental frequency of the nonlinear 
problem is not So but B = 8 o + e6 i+ . . . ; however, B must reduce 
to Bo when E$ 2 =0. This frequency shift can account for the 
phenomenon of "tracking" that has been observed experimentally 
[2] . Similarly, the frequency y = <J>-Bo that appears in the 
expression for |W 0 | should be considered as y = <J>-B. Then 
1 / 1 W o | shows all frequencies ny and W 0 /|W 0 | shows all fre- 
quencies ny±B, for n=0 ,1 , . . . This suggest that M has a com- 
plex Fourier series of the form 


l s n exp(inyTx). 
n=-°° 

Another factor of M must also be included since numerical 
examples show that M=0 if E<J> 2 is greater than some fixed 
value. This is similar tc the behavior of the van der Pol 
oscillator; see [6]. Thus, one may speculate that M has a 
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factor of the form F=l/[l+exp(-nTi) ] where n=n(E<|> 2 ). 
This would imply that F+l as when n^O and F-+0 as 
when n<0. Thus, M looks like 

oo 

1/ [l+exp(-nTi) 1 I s n exp (inyTi). 
n= -» 


PSD plots of R show frequencies of ny only while plots 
of Y show frequencies of <p and ny±B. 

3 . Examples : 

Example 1. In this example the system constants 
used are these: u=0, m=l lb.-s. 2 /in., C s =240 lb.-s./in., 

K s =0 . , Kb*l ,305,000 lb. /in., Q s = 200,000 lb. /in., 

6=. 0000285 in., and w=500 Hertz =100011 rad./s. Thus, 

Bo=833.33 rad./s. and a=. 000060915 in. The system is made 
nondimens ional using a for the g -displacement and B for the 
a-frequency. With these choices, the constants of this 
equation 

W’ ’ + CW '+[k(l-A/|w|)-iB]W=E4> 2 exp(i4>T) 

have these values: C=.288 , k=1.8792, A=. 467865 , B=.288 , and 
4> = 6tt / 5 . 

Figures 1 and 2 show changes in the solution Y vs . Z as E 
assumes the values 100n/ (IOOO tt) 2 a for n=0,l,...,7. The 
graphs are plotted for .2<t£.5s. The initial circle 
(for E=0) opens into an annular region, which becomes larger 
and thicker as E increases until a (transition) value of E 
occurs and the coefficient of exp (iBr) becomes zero. Thus, 
W=Nexp (i4»x ) , a circle of radius | N { . As E increases beyond 
this transition value, the solution remains a circle 
(figure 2.d) with radius | N | = | E4> 2 / ( -4> 2 + ic4> + k (1 - A/ j N| ) - iB. | . 

The characteristic, nonlinear frequency Bo is the angular 
frequency of the circle when E=0 , but this frequency increases 
as E increases. This tracking phenomenon is displayed in the 
PSD plots of figures 3 and 4 using the dimensional frequency B. 
In these figures only the 120-180 Hertz range is shown. At 
E=7/10 ,000II 2 a , there is no frequency in this range; instead, 
the circle is tranversed at the forcing frequency, u. 
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Figure 5, a typical full PSD plot, is the case 
H-4/10 ,000n 2 a. As shown earlier, one expects frequencies 
of u> and 6 to appear, as well as harmonics of ny ± 3 
where y*oj-B and n*l , 2, ... Thus, with 6*150 Hertz, and 
u)-500Hertz, one predicts that the PSD plot will exhibit 
peaks at 150, 200, 500, 550, 850, 900, ... Hertz. 

Example 2. In this example, the system constants 
from CDC[3] are these: y*0, m*. 20422 lb.-s. 2 /in., 

C s “20u lb.-s./in., K § *200,000 lb. /in., KuM , 000,000 lb. /in., 
oj* 500 Hertz, Q s =C s u)/2 lb. /in., and 6*. 0005 in. Thus, 

8*250 Hertz*500-rTrad . /s . and a*. 0007183 in. Figure 6 sum- 
marizes changes in Y vs. Z as E varies by 2S0n/ (IOOOtt) 2 a , 
n*0,l,...,4. These graphs are plotted for .45<t <.5s. 

In general, if the time interval is sufficiently large, 
the plots of annular regions will be completely filled 
(visually if not mathematically). However, when the ratio 
of 8 to y is p/q for small, positive integers p and q, then 
the curve actually falls on itself as time evolves and some 
attractive patterns appear. Figure 6.d is a case where 
8/y = 5/4, and the picture would look essentially the same 
whether shown for .45<t£.5 s. (actual) or for .45<t_<50 s. 

On the other hand, figure 6.c is more typical and would 
show a black annulus for .45<t£50 s. 

Example 3. This example uses the same data as example 1, 
but considers variations in the forcing frequency rather 
than E; i.e., E* . 7 and 4> varies. 

For each frequency 4> one expects to see a plot similar 
to one of those shown for example 1; thus, at each value 
of 0 either a circle or an annulus will describe Y vs. Z. 
Furthermore, based on continuity considerations, the 4>-axis 
will be subdivided into intervals within which each curve 
is a circle or an annulus. 

Example 1 showed that as E increased, the magnitude 
of the forcing function, B<p 2 , increased, and all plots 
were circles above the transition ^alue of E (and hence 
E<f> 2 ). A more careful analysis shows that this transition 
value is based not on the magnitude of the forcing function, 
but rather on the magnitude of the response at the forcing 
frequency; i.e., on E4> 2 / (-4» 2 + iC4> + <; • In example 1 , 4> is far 
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from Bo and the magnitude of the response is roughly E. 

But in the present example, one must consider values of 
<j> near and far from Bo. For those values near Bo, not 
only tie size of E but also the proximity of $ to Bo 
will make the magnitude of the response large. Therefor'?, 
one may predict that a circle at the forcing frequency will 
be seen for a plot of Y vs. Z when <P is near B, except for 
very small values of E. Elsewhere, one predicts an annular 
region unless the value of E is large. 

Figures 8 and 9 show the response values Y vs . Z for 
.4<tj<.5 s. as $ assumes the values .33, .385 , .66 , .99 , 

1.10, 1.32, 1.43, 1.65, 2.2, and 3.3 with E held at .7. 

Here the o-frequency for nondimensionalizing is w 0 . These 
irregularly-spaced frequencies were chosen for these reasons: 

a. .33 and .385 bound the value of transition from 

annulus to circle; 

b. .66, .99, 1.10 are values close to B; 

c. 1.32 and 1.43 bound the value of transition from 

circle to annulus; 

d. 1.65 has a nice picture; 

e. 2.20 and 3.30 are limiting values of interest 

At another value of E, the pictures would be similar although 
the transition values would be different. The extreme case 
E=0 is a circle of radius a with frequency B 0 , regardless of 
d> T s value. This was shown in figure l.a. 

If a three-dimensional plot Y vs. Z vs. 4>(=u>/u)o) were 
made for fixed E, one would find a frequency-response plot 
similar to the usual plots in linear theory. In particular, 
the maximum value of (Y z +Z 2 ) 1 ' 2 occurs near <#> * 1 . This is 
shown in figure 8.d where :he scale is four times as large 
as the scale of the other plots. The height of the response 
curves corresponds to the radius of the steady-state response 
when the plot is a circle or the outer radius when the plot 
is an annulus. In general, these response curves will be 
discontinuous at the transition points. 
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4. Conclusions: 


In studying the Jeffcott rotor with deadband and 
sinusoidal forcing, one must consider these three fre- 
quencies: (a) the forcing frequency w; (b) the natural 
frequency, u 0 » of the associated linear problem (deadband 
■6*0) ; and (c) the characteristic, nonlinear frequency Bo. 

The frequency u depends only on the forcing function; w 0 
depends only on the system parameters; B, with its base 
value Bo, depends on both the foiring function and the 
system parameters. 

For a given system and a nonzero, external, sinusoidal 
force, the y-z response is either a circle at the forcing 
frequency or an annulus composed of the (major) frequencies 
B and w as well as the (minor) harmonic frequencies 
n(u)-B)i8, for positive integers n. 

Frequency-response curves for a nonlinear problem and 
its associated (6“0) problem are similar, both reaching 
a maximum near w*u)o. Each point of the nonlinear plot, 
however, may represent either the radius of a circle or 
the outer radius of an annulus. There are also jump dis- 
continuities in the response curves at points of transition 
between circles and annuli. 

Two related problems, which show similar results in 
preliminary analyses, are rotors with sideforces and rotors 
with rubbing. Since the analysis of section 2 remains 
valid if E<J> ? is replaced by one coefficient G, one may 
consider sideforces as sinusoidals of the form GexpCiiox) 
where , A rubbing problem may be reduced tc a problem 
with a sideforce if one examines the response in the rotating 
rather than inertial coordinate system. Both of these pro- 
blems will be considered more carefully in future work. 

The principal obstacle remaining in this analysis is 
that of finding an explicit expression for transition values; 
i.e., expressing the transition point from an annulus to a 
circle as a function of the system and forcing parameters. 

This analysis does not allow two or more external forces 
(e.g., mass imbalance and sideforces). These problems appear 
to introduce no ’.ew theory, but do increase the computational 
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complications since one must consider not only forcing 
frequencies wi,a ) 2 , but also harmonic frequencies 

(1)X+0J2» Wl-(A>2, etc. 

Theaddition of asymmetric stiffness introduces diffi- 
culties that promise intriguing analyses based cn prelim- 
inary Runge-Kutta solutions. Not only do the circle/annulus 
plots become elliptic and occur with their axes rotated with 
respect to the y, z axes, but there may be other shapes and 
more than one transition point to consider. These problems, 
however, greatly extend the model's mimicry of an observed 
rotor's behavior. 

Finally, one needs to cons der stability questions for 
these nonlinear problems. Superficially, their stabilities 
are dictated by the corresponding linear problem's stability, 
but this has not yet been proven. 
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